1 Overview

This script:

  • Part 1: visualizes and analyzes the mean amplitude of response (meanbeta) in 8 focal regions (left and right SMG, left and right STS, left and right V1, bilateral APC and RFC).
  • Part 2: does the same across other, non-focal ROIs.
  • Part 3: visualizes and analyzes the reliability of domains with events and events within domains across all ROIs. This is the analog to the pre-registered MVPA analysis that is reported in the SI.

It also contains the code needed to reproduce Figure 3 in the main text of the paper. All instances of write_rds, write_csv, etc have been commented out to avoid writing over existing files. Uncomment those statements in order to reproduce these steps and save (new) outputs.

1.1 Read in data

univariate_data <- readRDS(here("outputs/univariate_data/study1_univariate_data_singlebetas.Rds"))
univariate_data$event <- relevel(as.factor(univariate_data$event), ref = "unexp")

univariate_summary_domain <- readRDS(here("outputs",
                                          "univariate_descriptive_summaries",
                                          "study1_univariate_summary_domain_allvoxelN.Rds"))

univariate_summary_task <- readRDS(here(
                                "outputs",
                            "univariate_descriptive_summaries",
 "study1_univariate_summary_task_allvoxelN.Rds"))

2 PART 1: Focal regions

focal_data <- univariate_data %>%
  filter(focal_region == 1)

focal_summary_domain <- univariate_summary_domain %>%
  filter(focal_region == 1)

focal_summary_task <- univariate_summary_task %>%
  filter(focal_region == 1)

domain_specific_regions <- c("physics", "psychology")
domain_general_regions <- c("early_visual", "MD")
focal_data_100 <- focal_data %>%
  filter(top_n_voxels == "100")

focal_summary_domain_100 <- focal_summary_domain %>%
  filter(top_n_voxels == "100")

focal_summary_task_100 <- focal_summary_task %>%
  filter(top_n_voxels == "100")

2.1 Does the VOE effect diminish in size over runs?

model_habituation <- lmer(data = focal_data_100 %>%
                            filter(ROI_category != "early_visual",
                                   event != "fam",
                                   ROI_category != "MD"), # cannot model MD regions because of non-independence of ROI selection (based on runs 2-4)
                          formula = meanbeta ~  extracted_run_number * event + (1|subjectID))
check_model(model_habituation)

plot(allEffects(model_habituation))

summary(model_habituation)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ extracted_run_number * event + (1 | subjectID)
##    Data: focal_data_100 %>% filter(ROI_category != "early_visual", event !=  
##     "fam", ROI_category != "MD")
## 
## REML criterion at convergence: 17809
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8142 -0.5931 -0.0313  0.5983  5.9668 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.2364   0.4862  
##  Residual              3.4448   1.8560  
## Number of obs: 4352, groups:  subjectID, 17
## 
## Fixed effects:
##                                Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)                   3.181e-01  1.212e-01  1.600e+01   2.624  0.01843
## extracted_run_number1         5.679e-02  4.873e-02  4.328e+03   1.165  0.24396
## extracted_run_number2        -9.594e-02  4.873e-02  4.328e+03  -1.969  0.04904
## extracted_run_number3         7.969e-03  4.873e-02  4.328e+03   0.164  0.87010
## event1                        7.583e-02  2.814e-02  4.328e+03   2.695  0.00706
## extracted_run_number1:event1  1.519e-01  4.873e-02  4.328e+03   3.117  0.00184
## extracted_run_number2:event1 -1.063e-02  4.873e-02  4.328e+03  -0.218  0.82738
## extracted_run_number3:event1 -6.623e-02  4.873e-02  4.328e+03  -1.359  0.17417
##                                
## (Intercept)                  * 
## extracted_run_number1          
## extracted_run_number2        * 
## extracted_run_number3          
## event1                       **
## extracted_run_number1:event1 **
## extracted_run_number2:event1   
## extracted_run_number3:event1   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) ext__1 ext__2 ext__3 event1 e__1:1 e__2:1
## extrctd_r_1  0.000                                          
## extrctd_r_2  0.000 -0.333                                   
## extrctd_r_3  0.000 -0.333 -0.333                            
## event1       0.000  0.000  0.000  0.000                     
## extrct__1:1  0.000  0.000  0.000  0.000  0.000              
## extrct__2:1  0.000  0.000  0.000  0.000  0.000 -0.333       
## extrct__3:1  0.000  0.000  0.000  0.000  0.000 -0.333 -0.333
model_habituation_results <- cbind(gen.m(model_habituation), gen.ci(model_habituation)[3:10,])
## Computing profile confidence intervals ...
model_habituation_results_byrun <- lsmeans(model_habituation, pairwise ~ event | extracted_run_number, pbkrtest.limit = 4352)$contrasts %>% as.data.frame()
DT::datatable(model_habituation_results, options = list(scrollX = TRUE))
DT::datatable(model_habituation_results_byrun, options = list(scrollX = TRUE))
# write_rds(model_habituation_results, path = here("outputs/univariate_results/habituation_over_runs_modelsummary.Rds"))
# write_rds(model_habituation_results_byrun, path = here("outputs/univariate_results/habituation_over_runs_perrun.Rds"))

2.2 Visualize by domain

regions <- levels(as.factor(focal_data$ROI_name_final))

plot_by_domain_run1 <- function(region, ymin, ymax) {
  plotobject <- ggplot(data = focal_summary_domain_100 %>% filter(str_detect(ROI_name_final, region), extracted_run_number == "run1", event != "fam"),
           aes(x = event, y = meanbeta, fill = domain)) +
      geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
      geom_errorbar(
        aes(ymin = meanbeta - se, ymax = meanbeta + se),
        position = position_dodge(width = .9),
        width = .2,
        colour = "black"
      ) +
      theme_cowplot(10) +
      facet_wrap( ~ ROI_name_final + domain, nrow = 1) +
      scale_fill_manual(values = c("#00AFBB", "#FC4E07")) +
      ylab("Amplitude") +
      xlab("Event") +
      theme(axis.text.x = element_text(
        angle = 90,
        vjust = 0.5,
        hjust = 1
    ),
    legend.position = "none") +
    coord_cartesian(ylim = c(ymin, ymax)) +
    ggtitle(paste0("ROI:", region))
  
  plotobject
}

This plot is a subset of Figure 3 in the paper.

SMG <- plot_by_domain_run1("supramarginal_", -1, 4) 

STS <- plot_by_domain_run1("superiortemporal_", 0, 4)

early_vis <- plot_by_domain_run1("V1_", 0, 11)

MD_APC <- plot_by_domain_run1("antParietal_bilateral", -1, 1)

MD_RFC <- plot_by_domain_run1("precentral_inferiorfrontal_R", -1, 1)

(SMG + STS) / (early_vis + MD_APC + MD_RFC) + plot_layout(widths = c(1, 1, 2, 1))
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

2.2.1 Plot run by run (except for MD regions)

regions <- levels(as.factor(focal_data$ROI_name_final))

plot_by_domain_all_runs <- function(region, ymin, ymax) {
  plotobject <- ggplot(data = focal_summary_domain_100 %>% filter(str_detect(ROI_name_final, region)),
           aes(x = event, y = meanbeta, fill = domain)) +
      geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
      geom_errorbar(
        aes(ymin = meanbeta - se, ymax = meanbeta + se),
        position = position_dodge(width = .9),
        width = .2,
        colour = "black"
      ) +
      theme_cowplot(10) +
      facet_wrap( ~ ROI_name_final + top_n_voxels + domain + extracted_run_number, nrow = 1) +
      scale_fill_manual(values = c("#00AFBB", "#FC4E07")) +
      ylab("Average beta (amplitude)") +
      xlab("Event") +
      theme(axis.text.x = element_text(
        angle = 90,
        vjust = 0.5,
        hjust = 1
    ),
    legend.position = "none") +
    coord_cartesian(ylim = c(ymin, ymax)) +
    ggtitle(paste0("ROI:", region))
  
  plotobject
}
plot_by_domain_all_runs("superiortemporal_L", -1, 3) | plot_by_domain_all_runs("superiortemporal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_by_domain_all_runs("supramarginal_L", -1, 1) | plot_by_domain_all_runs("supramarginal_R", -1, 1)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_by_domain_all_runs("V1_L", -1, 7) | plot_by_domain_all_runs("V1_R", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_by_domain_run1("antParietal_bilateral", -1, 1) | plot_by_domain_run1("precentral_inferiorfrontal_R", -1, 1)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

2.3 Visualize by task

regions <- levels(as.factor(focal_summary_task_100$ROI_name_final))

plot_univar_bytask <- function(region, ymin, ymax) {
  plotobject <- ggplot(data = focal_summary_task_100 %>% 
                         filter(str_detect(ROI_name_final, region),
                                extracted_run_number == "run1"),
           aes(x = event, y = meanbeta, fill = task)) +
      geom_bar(stat = "identity", aes(alpha = event), colour = "black") +
      geom_errorbar(
        aes(ymin = meanbeta - se, ymax = meanbeta + se),
        position = position_dodge(width = .9),
        width = .2,
        colour = "black"
      ) +
      # geom_point(data = focal_data_100 %>% filter(str_detect(ROI_name_final, region),
      #                                              extracted_run_number == "run1"),
      #            alpha = .1) +
      # geom_line(data = focal_data_100 %>% filter(str_detect(ROI_name_final, region)),
                # aes(group =
                      # subjectID),
                # alpha = .1) +
      theme_cowplot(10) +
      facet_wrap( ~ ROI_name_final + factor(task, c("solidity", "support", "goal", "efficiency")), nrow = 1) +
      scale_fill_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
      ylab("Average beta (amplitude)") +
      xlab("Event") +
      theme(axis.text.x = element_text(
        angle = 90,
        vjust = 0.5,
        hjust = 1
    ),
    legend.position = "none") +
    coord_cartesian(ylim = c(ymin, ymax)) +
    ggtitle(paste0("ROI:", region))
  
  plotobject
}
plot_univar_bytask("superiortemporal_L", -1, 3) | plot_univar_bytask("superiortemporal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_bytask("supramarginal_L", -1, 3) | plot_univar_bytask("supramarginal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_bytask("V1_L", -1, 7) | plot_univar_bytask("V1_R", -1, 7)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

plot_univar_bytask("antParietal_bilateral", -1, 3) | plot_univar_bytask("precentral_inferiorfrontal_R", -1, 3)
## Warning: Using alpha for a discrete variable is not advised.
## Using alpha for a discrete variable is not advised.

2.4 Results

2.4.1 Model fitting per domain

effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)

rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "domain1", "rep.L", "event:domain")

data_formodeling <- univariate_data %>%
  filter(event != "fam")

data_formodeling %>%
  group_by(top_n_voxels) %>%
  summarise(n = n())


ROIs <- levels(as.factor(data_formodeling$ROI_name_final))
top_n_voxels <- unique(data_formodeling$top_n_voxels)
runs <- c("run1", "run2", "run3", "run4", "all_runs")
modelsummaries_init <- data.frame()

for (ROI in ROIs) {
  curROI <- ROI
  for (topN in top_n_voxels) {
    cur_top_voxels_N <- topN
    
    for (run in runs) {
      curRun <- run
      
      if (curRun != "all_runs") {
        ROIdata <- data_formodeling %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN,
                 extracted_run_number == curRun)
        model <- lmer(data = ROIdata,
                      formula = meanbeta ~ event * domain + rep + (1 |
                                                                     subjectID))
        
      } else if (curRun == "all_runs") {
        ROIdata <- data_formodeling %>%
          filter(ROI_name_final == curROI,
                 top_n_voxels == topN)
        model <- lmer(data = ROIdata,
                      formula = meanbeta ~ event * domain + rep + (1 |
                                                                     subjectID))
        
      }
      
      model_just_maineffects <-
        update(model, formula = ~ . - event:domain)  # Without interaction term
      model_just_domain <-
        update(model_just_maineffects, formula = ~ . - event)
      model_just_event <-
        update(model_just_maineffects, formula = ~ . - domain)
      
      model_no_rep <- update(model, formula = ~ . - rep)
      
      
      BF_BIC_interaction <-
        data.frame(BF = exp((
          BIC(model_just_maineffects) - BIC(model)
        ) / 2),
        BF_interpretation = interpret_bf(exp((
          BIC(model_just_maineffects) - BIC(model)
        ) / 2)))  # BICs to Bayes factor
      
      BF_BIC_event <-
        data.frame(BF = exp((
          BIC(model_just_domain) - BIC(model_just_maineffects)
        ) / 2),
        BF_interpretation = interpret_bf(exp((
          BIC(model_just_domain) - BIC(model_just_maineffects)
        ) / 2)))
      
      BF_BIC_domain <-
        data.frame(BF = exp((
          BIC(model_just_event) - BIC(model_just_maineffects)
        ) / 2),
        BF_interpretation = interpret_bf(exp((
          BIC(model_just_event) - BIC(model_just_maineffects)
        ) / 2)))
      
      BF_BIC_rep <-
        data.frame(BF = exp((BIC(model_no_rep) - BIC(model)) / 2),
                   BF_interpretation = interpret_bf(exp((
                     BIC(model_no_rep) - BIC(model)
                   ) / 2)))
      
      BF_initial <- rbind(BF_BIC_event,
                          BF_BIC_domain,
                          BF_BIC_rep,
                          BF_BIC_interaction)
      
      rownames(BF_initial) <- rownames_no_intercept
      BFs <- rbind(BF_intercept_row,
                   BF_initial)
      
      effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
        select(d)
      effect_sizes <-
        rbind(effect_sizes_intercept_row, effect_sizes_initial)
      
      if (ROIdata$ROI_category[1] %in% domain_specific_regions) {
        ROI_category <- "domain_specific"
      } else {
        ROI_category <- "domain_general"
      }
      
      ROI_focal <- ROIdata$focal_region[1]
      modelsummary <-
        cbind(
          cur_top_voxels_N,
          curRun,
          ROI_category,
          ROI_focal,
          curROI,
          rownames(summary(model)$coefficients),
          summary(model)$coefficients,
          confint.merMod(model)[3:7, ],
          effect_sizes,
          BFs
        ) %>%
        as.data.frame()
      
      colnames(modelsummary) <-
        c(
          "top_n_voxels",
          "extracted_run_number",
          "domain",
          "focal_region",
          "ROI",
          "effect",
          "B",
          "SE",
          "df",
          "t",
          "p",
          "CI_95_lower",
          "CI_95_upper",
          "d",
          "BF",
          "BF_interpretation"
        )
      modelsummaries_init <-
        rbind(modelsummaries_init, modelsummary)
      
      cat(
        tab_model(
          model,
          show.std = TRUE,
          show.se = TRUE,
          show.stat = TRUE,
          auto.label = TRUE,
          title = paste0("Univariate ROI results in: ", curROI)
        )$knitr,
        "  \n  \n"
      )
    }
  }
}
modelsummaries <- modelsummaries_init %>%
  mutate_at(c(7:15), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

modelsummaries$effect <-
  as.factor(modelsummaries$effect)

levels(modelsummaries$effect) <-
  c("intercept", "domain", "event", "event:domain", "rep")
rownames(modelsummaries) <- NULL
# write_rds(modelsummaries, here("outputs", "univariate_results", "modelsummaries_allregions_alltopN.Rds"))
focal_modelsummaries <- read_rds(here("outputs/univariate_results/modelsummaries_allregions_alltopN.Rds")) %>%
  filter(focal_region == 1) 

DT::datatable(focal_modelsummaries %>% filter(top_n_voxels == "100", focal_region == 1) %>% arrange(effect), options = list(scrollX = TRUE,
                                                                                                                            pageLength = 100))
  # filter(extracted_run_number == "all_runs" | extracted_run_number == "run1") %>%

2.4.1.1 LSMG

LSMG_model <- lmer(data = univariate_data %>%
                     filter(top_n_voxels == "100",
                            event != "fam",
                            extracted_run_number == "run1", 
                            ROI_name_final == "supramarginal_L"),
                   formula = meanbeta ~ event * domain + rep + (1|subjectID))

LSMG_results_bydomain_top100 <- lsmeans(LSMG_model, pairwise ~ event | domain)$contrasts %>% as.data.frame()

DT::datatable(LSMG_results_bydomain_top100, options = list(scrollX = TRUE))
# write_rds(LSMG_results_bydomain_top100, path = here("outputs/univariate_results/LSMG_results_bydomain_top100.Rds"))

2.4.1.2 RSMG

RSMG_model <- lmer(data = univariate_data %>%
                     filter(top_n_voxels == "100",
                            event != "fam",
                            extracted_run_number == "run1", 
                            ROI_name_final == "supramarginal_R"),
                   formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(RSMG_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
##    Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",  
##     extracted_run_number == "run1", ROI_name_final == "supramarginal_R")
## 
## REML criterion at convergence: 1047.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2569 -0.5727  0.0203  0.6278  2.5903 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 1.517    1.232   
##  Residual              2.312    1.520   
## Number of obs: 272, groups:  subjectID, 17
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)     -0.01345    0.31266  16.00000  -0.043  0.96621   
## event1           0.21889    0.09219 251.00000   2.374  0.01834 * 
## domain1          0.28691    0.09219 251.00000   3.112  0.00207 **
## rep.L           -0.36357    0.13038 251.00000  -2.789  0.00570 **
## event1:domain1   0.05727    0.09219 251.00000   0.621  0.53505   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 rep.L
## event1      0.000                     
## domain1     0.000  0.000              
## rep.L       0.000  0.000  0.000       
## event1:dmn1 0.000  0.000  0.000  0.000
lsmeans(RSMG_model, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.552 0.261 251 2.118   0.0352 
## 
## domain = psychology:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.323 0.261 251 1.240   0.2163 
## 
## Results are averaged over the levels of: rep 
## Degrees-of-freedom method: kenward-roger

2.4.1.3 LSTS

LSTS_model <- lmer(data = univariate_data %>%
                     filter(top_n_voxels == "100",
                            event != "fam",
                            extracted_run_number == "run1", 
                            ROI_name_final == "superiortemporal_L"),
                   formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(LSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
##    Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",  
##     extracted_run_number == "run1", ROI_name_final == "superiortemporal_L")
## 
## REML criterion at convergence: 1051.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4985 -0.5627 -0.0406  0.6461  2.9036 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 1.773    1.331   
##  Residual              2.327    1.525   
## Number of obs: 272, groups:  subjectID, 17
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      0.54190    0.33591  16.00000   1.613 0.126244    
## event1           0.16841    0.09249 251.00000   1.821 0.069839 .  
## domain1         -0.32878    0.09249 251.00000  -3.555 0.000452 ***
## rep.L            0.13244    0.13081 251.00000   1.012 0.312286    
## event1:domain1   0.04951    0.09249 251.00000   0.535 0.592897    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 rep.L
## event1      0.000                     
## domain1     0.000  0.000              
## rep.L       0.000  0.000  0.000       
## event1:dmn1 0.000  0.000  0.000  0.000
lsmeans(LSTS_model, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.436 0.262 251 1.666   0.0970 
## 
## domain = psychology:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.238 0.262 251 0.909   0.3643 
## 
## Results are averaged over the levels of: rep 
## Degrees-of-freedom method: kenward-roger

2.4.1.4 RSTS

RSTS_model <- lmer(data = univariate_data %>%
                     filter(top_n_voxels == "100",
                            event != "fam",
                            extracted_run_number == "run1", 
                            ROI_name_final == "superiortemporal_R"),
                   formula = meanbeta ~ event * domain + rep + (1|subjectID))
summary(RSTS_model)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: meanbeta ~ event * domain + rep + (1 | subjectID)
##    Data: univariate_data %>% filter(top_n_voxels == "100", event != "fam",  
##     extracted_run_number == "run1", ROI_name_final == "superiortemporal_R")
## 
## REML criterion at convergence: 1015.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.16150 -0.59428 -0.03496  0.59256  2.87577 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 2.143    1.464   
##  Residual              1.992    1.411   
## Number of obs: 272, groups:  subjectID, 17
## 
## Fixed effects:
##                 Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)      0.95600    0.36522  16.00000   2.618 0.018662 *  
## event1           0.25840    0.08557 251.00000   3.020 0.002792 ** 
## domain1         -0.30885    0.08557 251.00000  -3.609 0.000371 ***
## rep.L           -0.22514    0.12102 251.00000  -1.860 0.064000 .  
## event1:domain1  -0.05083    0.08557 251.00000  -0.594 0.553030    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) event1 doman1 rep.L
## event1      0.000                     
## domain1     0.000  0.000              
## rep.L       0.000  0.000  0.000       
## event1:dmn1 0.000  0.000  0.000  0.000
lsmeans(RSTS_model, pairwise ~ event | domain)$contrasts
## domain = physics:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.415 0.242 251 1.715   0.0876 
## 
## domain = psychology:
##  contrast    estimate    SE  df t.ratio p.value
##  unexp - exp    0.618 0.242 251 2.555   0.0112 
## 
## Results are averaged over the levels of: rep 
## Degrees-of-freedom method: kenward-roger

2.4.2 Model fitting per task

effect_sizes_intercept_row <- data.frame(d = NA)
BF_intercept_row <- data.frame(BF = NA, BF_interpretation = NA)

rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames(BF_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1", "rep.L")

data_formodeling <- univariate_data %>%
  filter(event != "fam",
         top_n_voxels == 100)

data_formodeling %>%
  group_by(top_n_voxels) %>%
  summarise(n = n())

ROIs <- levels(as.factor(data_formodeling$ROI_name_final))
top_n_voxels <- unique(data_formodeling$top_n_voxels)
runs <- c("run1", "all_runs")

tasks <- unique(data_formodeling$task)
ROIs <- unique(data_formodeling$ROI_name_final)
data_formodeling$event <- relevel(data_formodeling$event, ref = "unexp")

modelsummaries_pertask_init <- data.frame()



for (ROI in ROIs) {
  curROI <- ROI
  for (task in tasks) {
    cur_task <- task
    
    for (run in runs) {
      if (run == "run1") {
        ROIdata <-
          data_formodeling %>% filter(ROI_name_final == curROI,
                                      task == cur_task,
                                      extracted_run_number == "run1")
      } else if (run == "all_runs")
        ROIdata <-
          data_formodeling %>% filter(ROI_name_final == curROI,
                                      task == cur_task)
      
        model <- lmer(data = ROIdata,
                      formula = meanbeta ~ event + rep + (1 | subjectID)) # had to remove run random int (convergence))

      focal_region <- ROIdata$focal_region[1]
      task_domain <- as.character(ROIdata$domain[1])
      ROI_category <- ROIdata$ROI_category[1]

      
      model_null <-
        update(model, formula = ~ . - event)  # Without interaction term
      
      model_norep <-
        update(model, formula = ~ . - rep)  # Without interaction term
  
      BF_BIC_event <-
        data.frame(BF = exp((BIC(model_null) - BIC(model)) / 2),
                   BF_interpretation = interpret_bf(exp((
                     BIC(model_null) - BIC(model)
                   ) / 2)))  # BICs to Bayes factor
      
      BF_BIC_rep <-
        data.frame(BF = exp((BIC(model_norep) - BIC(model)) / 2),
                   BF_interpretation = interpret_bf(exp((
                     BIC(model_norep) - BIC(model)
                   ) / 2)))  # BICs to Bayes factor      
      # BF_initial <- rbind(BF_BIC_event)
      
      
      BFs <- rbind(BF_intercept_row,
                   BF_BIC_event,
                   BF_BIC_rep)
      
      rownames(BFs)[2:3] <- rownames_no_intercept
      
      effect_sizes_initial <- lme.dscore(model, ROIdata, "lme4") %>%
        select(d)
      effect_sizes <-
        rbind(effect_sizes_intercept_row, effect_sizes_initial)
      cur_top_voxels_N <- 100

      modelsummary <-
        cbind(
          cur_top_voxels_N,
          run,
          ROI_category,
          focal_region,
          curROI,
          cur_task,
          rownames(summary(model)$coefficients),
          summary(model)$coefficients,
          confint.merMod(model)[3:5, ],
          effect_sizes,
          BFs
        ) %>%
        as.data.frame()
      
      colnames(modelsummary) <-
        c(
          "top_n_voxels",
          "extracted_run_number",
          "domain",
          "focal_region",
          "ROI",
          "task",
          "effect",
          "B",
          "SE",
          "df",
          "t",
          "p",
          "CI_95_lower",
          "CI_95_upper",
          "d",
          "BF",
          "BF_interpretation"
        )
      
      
      modelsummaries_pertask_init <-
        rbind(modelsummaries_pertask_init, modelsummary)
      # model_parameters(model) %>% print_md()
      cat(
        tab_model(
          model,
          show.std = TRUE,
          show.se = TRUE,
          show.stat = TRUE,
          auto.label = TRUE,
          title = paste0(
            "Exploratory Univariate ROI results for the task ",
            cur_task,
            "in the region: ",
            curROI,
            "in ",
            curRun
          )
        )$knitr,
        "  \n  \n"
      )
    }
  }
}
modelsummaries_pertask <- modelsummaries_pertask_init %>%
  mutate_at(c(8:16), as.numeric) %>%
  mutate(star = as.factor(
           case_when(p < .001 ~ "***",
                     p < .01 ~ "**",
                     p < .05 ~ "*",
                     p < .1 ~ "~",
                     TRUE ~ " ")
         )) %>%
  mutate(p = round(p, digits = 3))

modelsummaries_pertask$effect <-
  as.factor(modelsummaries_pertask$effect)

levels(modelsummaries_pertask$effect) <-
  c("intercept", "event", "rep")
rownames(modelsummaries_pertask) <- NULL
# write_rds(modelsummaries_pertask, here("outputs", "univariate_results", "modelsummaries_allregions_pertask.Rds"))
modelsummaries_pertask <- read_rds(here("outputs", "univariate_results", "modelsummaries_allregions_pertask.Rds"))

2.4.2.1 Plot model outputs

focal_modelsummaries_pertask_forplotting <- modelsummaries_pertask %>%
  filter(top_n_voxels == "100",
         effect == "event",
         extracted_run_number == "run1",
         focal_region == 1)

focal_modelsummaries_pertask_forplotting$ROI <- factor(focal_modelsummaries_pertask_forplotting$ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "antParietal_bilateral", "precentral_inferiorfrontal_R", "V1_L", "V1_R"))


pertask_b <- ggplot(focal_modelsummaries_pertask_forplotting,
       aes(ROI, B, colour = task, group = task)) +
  geom_point() +
  geom_errorbar(
    aes(ymin = B - SE, ymax = B + SE, 
    colour = task),
    position = position_dodge(width = .9),
    width = .2
  ) +  
  # geom_line() +
  geom_hline(yintercept = 0) +
  facet_wrap( ~ task, nrow = 1) +
  scale_colour_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
  ylab("B and SE") +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  )) +
  coord_cartesian(ylim = c(-.5, .7))

perregion_b <- ggplot(focal_modelsummaries_pertask_forplotting,
       aes(task, B, colour = task, group = task)) +
  geom_point() +
  geom_errorbar(
    aes(ymin = B - SE, ymax = B + SE, 
    colour = task),
    position = position_dodge(width = .9),
    width = .2
  ) +  
  # geom_line() +
  geom_hline(yintercept = 0) +
  facet_wrap( ~ ROI, nrow = 1) +
  scale_colour_manual(values = c("#deaf1f", "#f34b00", "#00aa8b", "#4f8e00")) +
  ylab("B and SE") +
  theme(axis.text.x = element_text(
    angle = 90,
    vjust = 0.5,
    hjust = 1
  )) +
  coord_cartesian(ylim = c(-.5, .7))
pertask_b

perregion_b

2.5 Sensitivity of focal region results to top voxel criteria

focal_modelsummaries$top_n_voxels <- factor(focal_modelsummaries$top_n_voxels, levels = c("10", "50", "100", "150", "200", "250", "300"))
focal_modelsummaries$extracted_run_number <- factor(focal_modelsummaries$extracted_run_number, levels = c("run1", "run2", "run3", "run4", "all_runs"))


ggplot(focal_modelsummaries %>%
         filter(extracted_run_number != "all_runs",
                ROI != "antParietal_bilateral",
                ROI != "precentral_inferiorfrontal_R",
                effect != "intercept",
                effect != "rep"),
       aes(extracted_run_number, d, colour = top_n_voxels)) +
  geom_point() +
  geom_line(aes(group = top_n_voxels)) +
  ylab("Cohen's D") +
  facet_wrap(~ effect + factor(ROI, levels = c("supramarginal_L", "supramarginal_R", "superiortemporal_L", "superiortemporal_R", "V1_L", "V1_R")), nrow = 3) +
  # coord_cartesian(ylim = c(-.75, 0.75)) +
  geom_hline(yintercept = 0) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

3 PART 2: Other regions

nonfocal_modelsummaries <- read_rds(here("outputs/univariate_results/modelsummaries_allregions_alltopN.Rds")) %>%
  filter(focal_region == 0) 

DT::datatable(nonfocal_modelsummaries %>% 
                filter(top_n_voxels == "100") %>% 
                arrange(effect), 
              options = list(scrollX = TRUE,
                             pageLength = 100))

4 PART 3: Region by region univar effect size

Here we repeat the MVPA manyregions analysis, whose results are difficult to interpret given the lack of reliable multivariate signal towards unexpected events, despite clear univariate effects.

data_for_modeling <- univariate_data %>%
  filter(event != "fam",
         top_n_voxels == "100") %>%
  filter(extracted_run_number == "run1") %>%
  filter(!str_detect(ROI_name_final, "bilateral"))

ROIs <- levels(as.factor(data_for_modeling$ROI_name_final))

effect_sizes_intercept_row <- data.frame(d = NA)

rownames(effect_sizes_intercept_row) <- "(Intercept)"
rownames_no_intercept <- c("event1")
modelsummaries_init <- data.frame()

for (ROI in ROIs) {
  curROI <- ROI

  # FIRST COMPUTE EVENT EFFECT SIZES, FOR EACH DOMAIN
  physics_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           domain == "physics")
  
  psychology_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           domain == "psychology")
      
  # had to remove random intercept for runs (convergence))
  physics_model <- lmer(data = physics_data,
                formula = meanbeta ~ event + iqr_contrast + mean_luminance + mean_motion + mean_hsf + mean_lsf + mean_rect + mean_curv + (1|subjectID))
  
  psychology_model <- lmer(data = psychology_data,
                formula = meanbeta ~ event + (1|subjectID))
  
  physics_event_effect_size <- lme.dscore(physics_model, physics_data, "lme4") %>%
    select(d)
  psychology_event_effect_size <- lme.dscore(psychology_model, psychology_data, "lme4") %>%
    select(d)
  
  # SECOND COMPUTE DOMAIN EFFECT SIZES, FOR EACH EVENT
  expected_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           event == "exp")
  
  unexpected_data <- data_for_modeling %>%
    filter(ROI_name_final == curROI,
           event == "unexp")
  
  # had to remove random intercept for runs (convergence))
  exp_model <- lmer(data = expected_data,
                        formula = meanbeta ~ domain + (1 | subjectID))
  
  unexp_model <- lmer(data = unexpected_data,
                           formula = meanbeta ~ domain + (1 | subjectID))
  
  exp_domain_effect_size <-
    lme.dscore(exp_model, expected_data, "lme4") %>%
    select(d)
  unexp_domain_effect_size <-
    lme.dscore(unexp_model, unexpected_data, "lme4") %>%
    select(d)
  
  if (physics_data$ROI_category[1] == "MD" | physics_data$ROI_category[1] == "early_visual") {
    ROI_domain <- "general"
  } else {
    ROI_domain <- "specific"
  }
  
  modelsummary <-
    cbind(
      ROI_domain,
      physics_data$ROI_category[1],
      curROI,
      physics_event_effect_size,
      psychology_event_effect_size,
      exp_domain_effect_size,
      unexp_domain_effect_size
    ) %>%
    as.data.frame()
    
      modelsummaries_init <-
        rbind(modelsummaries_init, modelsummary)
}
colnames(modelsummaries_init) <-
  c("ROI_domain",
    "ROI_category",
    "ROI",
    "physics_event_effect_size",
    "psychology_event_effect_size",
    "exp_domain_effect_size",
    "unexp_domain_effect_size"
  )
rownames(modelsummaries_init) <- NULL


# write_rds(modelsummaries_init, here("outputs/univariate_results/univar_manyregions_results.Rds"))
univar_manyregions_results <- read_rds( here("outputs/univariate_results/univar_manyregions_results.Rds"))
plot1 <- ggplot(univar_manyregions_results, aes(psychology_event_effect_size,       physics_event_effect_size)) +
    geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  coord_cartesian(xlim = c(-0.5, 1), ylim = c(-0.5, 1)) +
  geom_point(aes(colour = ROI_category), size = 3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ROI_domain, scales = "free") +
  # coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
  scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
  xlab("VOE Effect Size (d), Physics Events") +
  ylab("VOE Effect Size (d), Psychology Events") +
  ggtitle(label = "Exp 1 VOE effects for each domain across regions") +
  geom_text_repel(aes(label = ROI), size = 3) +
    theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        legend.position = "none")


plot2 <- ggplot(univar_manyregions_results, aes(exp_domain_effect_size, unexp_domain_effect_size)) +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  coord_cartesian(xlim = c(-1.2, 1.2), ylim = c(-1.2, 1.2)) +
  geom_point(aes(colour = ROI_category), size = 3) +
  geom_smooth(method = "lm") +
  facet_wrap(~ROI_domain, scales = "free") +
  # coord_cartesian(xlim = c(-1.5,1.5), ylim = c(-1.5,1.5)) +
  scale_colour_manual(values = c("#fb00d3", "#f8bf00", "#4894ff", "#f71d00")) +
  xlab("Domain Effect Size (d), Expected Events") +
  ylab("Domain Effect Size (d), Unexpected Events") +
  ggtitle(label = "Exp 1 Domain effects for each event across regions") +
  geom_text_repel(aes(label = ROI), size = 3) +
    theme(axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14),
        legend.position = "none")
plot1
## `geom_smooth()` using formula = 'y ~ x'

plot2
## `geom_smooth()` using formula = 'y ~ x'

# instead of testing whether the linear relationship between x and y is 0, 
# we test for independence instead. H0 is that x and y are independent; F_{XY}(x,y) = F_X(x) F_Y(y).
observed_cors_ind <- univar_manyregions_results %>%
  group_by(ROI_domain) %>%
  summarise(
    cor_domain_within_event =
      np.cor.test(
        exp_domain_effect_size,
        unexp_domain_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$estimate,
    
    p_domain_within_event =
      np.cor.test(
        exp_domain_effect_size,
        unexp_domain_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$p.value,
    
    cor_event_within_domain =
      np.cor.test(
        physics_event_effect_size,
        psychology_event_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$estimate,
    
    p_event_within_domain =
      np.cor.test(
        physics_event_effect_size,
        psychology_event_effect_size,
        alternative = "two.sided",
        independent = TRUE
      )$p.value
  )
observed_cors_ind
# write_rds(observed_cors_ind, path = here("outputs/univariate_results/study1_univar_manyregions_observed_cor.Rds"))

4.0.1 Bootstrap stats

set.seed(2020)
## First define a function to work out the difference of corrs:
diff_corr <- function(data, indices) {
  data <- data[indices,]
  cor1 <-
    np.cor.test(DS_results$exp_domain_effect_size,
                DS_results$unexp_domain_effect_size,
                alternative = "two.sided",
                independent = TRUE,
                parallel = TRUE)$estimate
  cor2 <-
    np.cor.test(
      data$psychology_event_effect_size,
      data$physics_event_effect_size,
      alternative = "two.sided",
      independent = TRUE,
      parallel = TRUE)$estimate

    return(cor1 - cor2)
}
DS_results <- univar_manyregions_results %>%
  filter(ROI_domain == "specific")

DG_results <- univar_manyregions_results %>%
  filter(ROI_domain == "general")

4.0.2 Domain specific

doMC::registerDoMC(cores = 2)  # for running in parallel

# Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DS <- boot(
#   data = DS_results,
#   R = 4000,
#   statistic = diff_corr,
#   stype = "i"
# )

# saveRDS(res_boot_DS, here("outputs/univariate_results/study1_univariate_dsregions_4000perms_confirmatory.Rds"))

res_boot_DS <- read_rds(here("outputs/univariate_results/study1_univariate_dsregions_4000perms_confirmatory.Rds"))

## Retrieve the empirical 95% confidence interval:
ds_results <- append(boot.ci(res_boot_DS, type = "perc", conf = 0.95),boot.pval(res_boot_DS))

# saveRDS(ds_results, here("outputs/univariate_results/study1_univariate_dsregions_summary.Rds"))

4.0.3 Domain general

## Then apply a bootstrap procedure with 4000 draws (uncomment to reproduce):
# res_boot_DG <- boot(data = DG_results,
#                  R = 4000,
#                  statistic = diff_corr,
#                  stype = "i")

# saveRDS(res_boot_DG, here("outputs/univariate_results/study1_univariate_dgregions_4000perms_confirmatory.Rds"))
res_boot_DG <- readRDS(here("outputs/univariate_results/study1_univariate_dgregions_4000perms_confirmatory.Rds"))

plot(res_boot_DG)

## Retrieve the empirical 95% confidence interval:

dg_results <- append(boot.ci(res_boot_DG, type = "perc", conf = 0.95),boot.pval(res_boot_DG))
# saveRDS(dg_results, here("outputs/univariate_results/study1_univariate_dgregions_summary.Rds"))